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We analyze the electronic properties of bilayer graphene with Bernal stacking and a low concentration of 
adatoms. Assuming that the host bilayer lies on top of a substrate, we consider the case where impurities are 
adsorbed only on the upper layer. We describe non-magnetic impurities as a single orbital hybridized with car¬ 
bon’s p;. states. The effect of impurity doping on the local density of states with and without a gated electric field 
perpendicular to the layers is analyzed. We look for Anderson localization in the different regimes and estimate 
the localization length. In the biased system, the field induced gap is partially filled by strongly localized im¬ 
purity states. Interestingly, the structure, distribution and localization length of these states depend on the field 
polarization. 
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I. INTRODUCTION 

Graphene, in all its allotropic forms, is a material with ex¬ 
ceptional mechanical, electronic and thermal properties. Its 
discovery led to one of the most active fields in material sci¬ 
ence and condensed matter research during the last decade. 
Graphene monolayer, usually referred simply as graphene, 
and multilayers have different properties due to a subtle dif¬ 
ference in their band structure. It is now well established 
that in graphene monolayers, the electronic excitations with 
crystal momentum close to the K or K' points of the Bril- 
louin zone (BZ), are chiral quasiparticles behaving as mass¬ 
less Dirac fermions. These excitations dominate the low tem¬ 
perature physics leading to a number of remarkable phenom¬ 
ena in clean samples.Impurities, adatoms and structural 
defects change these properties and there has been a consid¬ 
erable effort to study and characterize the different types of 
defects and disorder in graphene"*"*^ as well as the effect of 
doping on them. 

The problem of disorder and electron localization has at¬ 
tracted the attention of many groups for Dirac fermions 
tend to elude localization in systems with Anderson-type 
disorder. Impurities leading to short range disorder at the 
atomic scale generate inter-valley mixing and break the sym- 
plectic symmetry opening the route to strong localization. 

Bilayer graphene (BLG) presents some fundamental dif¬ 
ferences due to its crystallographic structure. It consist of 
a stacking of two graphene layers and in the most common 
structure, known as the Bernal stacking, only one of the two 
non-equivalent sites (A, E) of the honeycomb lattice of the top 
layer lies on top of a site of the bottom layer. The resulting 
structure, shown in Fig. 1, induces a weak coupling of the 
two layers. The unit cell has four carbon atoms leading to 
four TT-bands, two of them having a parabolic dispersion rela¬ 
tion around the K and K’ points of the BZ and that touch each 
other at the Fermi energy.^^^^^ 

In most of the experimental setups, BLG lies on top of a 
substrate and the impurities are adsorbed on the top layer only. 
When atoms like hydrogen or fluorine are adsorbed they are 
bounded to a single C atom. 

One of the most interesting aspects of this system is that 


its electronic structure can be controlled with an electric field 
applied perpendicular to the layers.In biased BLG a gap 
opens at the Fermi level and impurities may induce a bound 
state in the gap.^^’^^ As noted in Ref. [33], the impurity spec¬ 
tral density and the existence of the bound states may depend 
on the polarity of the field. A finite impurity concentration 
generates a gate dependent impurity band creating new and 
encouraging alternatives to control the transport properties. 
However, in contrast to the important activity in the study of 
disordered graphene, the problem of BLG with a diluted con¬ 
centration of adatoms inducing short range potentials has not 
been investigated in detail. 

In this work we study the problem of a low concentration of 
impurities in biased and unbiased BLG. We present a model 
that aims to describe fluorinated BLG, an extension of the 
model of Ref. [24] used to discuss the experiments of Ref. 

[13]. 

In section II we present the model and revisit the single im¬ 
purity problem. In section III we describe the numerical meth¬ 
ods and present results for the local density of states (LDOS) 
at the different sites. Section IV includes a discussion of lo¬ 
calization and a summary and conclusions are presented in 
section V. 


II. THE MODEL 

The Hamiltonian of the system is H - H^lg + //imp + //hyb 
where the first term describes the electronic structure of the 
BLG, the second one is the impurities’ Hamiltonian and the 
last one includes the hybridization between each impurity or¬ 
bital and the p^ orbital of the underlying C atom. In the tight- 
binding approximation the BLG Hamiltonian reads 

HbLC = - X 

i,k,cr 

^ikcT^ikcT /fccr) 

i,k,cr 

k,cr 
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FIG. 1. (a) Bernal structure of bilayer graphene. Panels (b) and (c) 
show the low energy band structure for the unbiased and biased BLG, 
respectively. In (d) the LDOS for the unbiased case, continuous line 
corresponds to the Bi and sites, dashed line to the A\ and 62 sites, 
(e) LDOS for a positive (V > 0) bias voltage. Continuous (dashed) 
line corresponds to the Bi (Ai) sublattice. For a negative voltage the 
LDOS can be obtained from the ones with positive V by replacing ai 
by -o). 
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FIG. 2. Color maps (in logarithmic scale) of the impurity spectral 
densities in the [oj, V] plane, (a) and (b) correspond to one impu¬ 
rity on the Ai and Bi sublattice, respectively. Dark areas indicate the 
electric field induced gap in the pristine sample. Lower panels show 
the impurity spectral density for the unbiased (dashed lines) and bi¬ 
ased (continuous lines) cases: left and right columns correspond to 
impurities on the Ai and B[ sublattices, respectively, and different 
polarities (as indicated in the insets) 


Here, and ( 7 ,^ 0 - destroy electrons with wavevector k and 
spin cr in sub-lattices A and B, respectively, and the subindex 
i - 1(2) refers to the top (bottom) plane. V is the bias voltage, 
t and G are the intra-plane and inter-plane hoppings, respec¬ 
tively, and (p{k) - YjS ^ where {^) are the three vectors 
connecting one site with its neighbors in the same plane. In 
our notation, the C atoms of the top layer in the A sublattice, 
referred as the A i sublattice, lie on top of the C atoms in the 
B sublattice of the bottom layer {B 2 sublattice), see Fig. 1. 

We describe non-magnetic impurities as single orbital im¬ 
purities where the electron-electron interactions are not ex¬ 
plicitly included, 

^imp = ^ > (2) 

1 , 0 - 

where creates an electron on the impurity orbital at site I 
and energy eo and the sum runs over the sites of carbon lat¬ 
tice having an adsorbed impurity on top. The last term of the 
Hamiltonian describes the hybridization of the impurity and 
the graphene orbitals of the top layer 

^^hyb = r 2 + 4iJio-) + r X ’ 

IgA.\,o’ /c5[,o” 

(3) 

with the sum taken over all the sites with an impurity on top, 
where R/ is the coordinate of site 


/—a similar expression holds for b^i^. Typical values of the 
microscopic parameters are f = 2.8 eV and = O.lf, while 
the bias voltage V is taken in the range |y| < 0.3 eV and with 
any loss of generality we take sq > 0. With our one electron 
Hamiltonian, the case eo < 0 can be obtained from the previ¬ 
ous one by an electron-hole transformation, i.e. by replacing a> 
by -a> and V by -V. In what follows we take eo = 0.3 eV and 
y = 5.6 eV. Figure 1 illustrates the BLG lattice and its band 
structure. As we will not consider spin dependent effects, we 
drop the spin index in what follows. 

It is instructive to review some aspects of the single impu¬ 
rity problem before presenting the many impurities case.^^’^^ 
For one impurity, the retarded impurity propagator = 
((/;>//)) takes the form 


to + i0+ - eo - '^{to, V) ’ 

where S((u, V) - y^gCto, V) is the impurity’s self-energy and 
g(a>, V) is the local propagator of electrons in the C orbital 
hybridized with the impurity. The carbon-carbon propagator 
can be evaluated in the continuous limit. Color maps of the 
impurity spectral densities for impurities on Ai and on Bi sites 
are shown in Fig. 2. 

Let us consider first the case of an impurity adsorbed on 
top of an Ai site (left panels of Fig. 2). For V = 0 the im- 
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purity spectral density pimp(w) = -l/7rlm^/; has the charac¬ 
teristic structure of a resonant state, the real part of the self¬ 
energy shifts the maximum from eo ^ 0 towards the Dirac 
point generating a narrow resonance at the renormalized en¬ 
ergy eo = eo + ReS(eo, 0). For the impurity parameters used in 
this calculation, the renormalized energy is an order of mag¬ 
nitude smaller than the bare energy (eq ^ sq). Close to the 
Dirac point the impurity spectral density shows the character¬ 
istic |w| behavior of the LDOS of the Ai sites. For V 0 a 
gap opens at the Dirac point and the impurity spectral density 
Pinjp(w) may show a bound state within the gap. An important 
effect of the polarity of the field V is apparent from the figure: 
a negative voltage V leads to a bound state within the energy 
gap close to the top of the valence band while for a positive V 
the bound state energy—if observed—lies exponentially close 
to the conduction band edge. This is due to the structure of 
the LDOS at the Ai sublattice: for positive V the LDOS at 
the edge of the conduction band Ec behaves as w - Ec, as in 
the 4D electron gas where a strong coupling to the impurity is 
required to split a bound state out of the band (see Fig. le)). 

For impurities on the Bi sites (right panels of Fig. 2) the 
results are somewhat different. For V = 0 the width of the 
impurity resonance is much broader due to the larger LDOS 
of the underlying C atom. For small and positive V the bound 
state lies close to the gap centre while for small negative V no 
bound state occurs. This effect can be understand by looking 
at the LDOS at the Bi sites in the biased BLG (see Fig le)). 
There, the LDOS of the conduction band for small and posi¬ 
tive V presents a ID-like van Hove singularity leading always 
to a bound state, while for V < 0 a 3D-like singularity at the 
edge of the conduction band requires a minimum value of the 
parameters to split a state out of the band. However, this ef¬ 
fect, discussed in Ref. [33], is observed only for extremely 
small values of the bias voltage. For physically relevant val¬ 
ues of the gap, bound states occur for both polarities although 
their position depends on the sign of V. 

These asymmetries illustrate the importance of the polarity 
of the electric field on the electronic structure of the impurity 
doped system. The variation of the impurity energy eo with V 
depends on the way the electric field is induced in the system 
and on the characteristic of the impurity. To minimize the 
number of parameters in the model we present results with 
constant, V independent, eo. Having in mind the one impurity 
problem results, summarized in Fig. 2, the more relevant case 
of many impurities can be easily interpreted. 


III. NUMERICAL RESULTS FOR THE MANY 
IMPURITIES CASE 


Here we present results for the case of a small concentra¬ 
tion of adatoms on the top layer. Calculations using Density 
Functional Theory show that, in the case of fluorine atoms, 
the adsorption energies on the two non-equivalent sites, Ai 
and Bi, are almost equal. Some estimations, however, sug¬ 
gest that there could be a tiny energy gain for adatoms on 
the Bi sites.Interestingly, in other carbon based systems 
like monolayer graphene with substitutional nitrogen impu- 
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FIG. 3. Impurity spectral densities for different bias voltages and 
different impurity distributions. Top [(a), (b), (c)], central [(d), (e), 

(f) ] and bottom [(g), (h), (i)] panels correspond to 100%, 75% and 
50% of the impurities on the Bi sublattice respectively. Left [(a), (d), 

(g) ], central [(b), (e), (h)] and right [(c), (f), (i)] columns correspond 
to L = 0, k = 0.1 eV and V = -O.leV, respectively. Vertical dashed 
lines indicate the gap corresponding to the pristine BLG. 


rities, it has been observed that as the nitrogen concentra¬ 
tion increases, impurities tend to be absorbed preferentially in 
one of the two equivalent sublattices.These self-organized 
structures of the nitrogen doped graphene are stabilized by 
the impurity-impurity interaction that favors impurities on the 
same sublattice, an effect that scales quickly with the impu¬ 
rity concentration.^® For the case of diluted fluorine adatoms 
on BLG there are no evidences of clustering on one sublattice. 
Moreover, the interaction between impurities on graphene is 
known to depend crucially on the type of impurity and on the 
adsorption geometry.Based on these facts, in what follows 
we consider different impurity distributions, going from 50% 
of the impurities in each sublattice to 100% of them on the 
Bi sites. We start with a detailed analysis of the LDOS at the 
impurity and at the different sites of the BLG. 


A. Spectral densities 

To calculate pimp(<y) in the many impurities case we 
first use the Chebyshev polynomials method which has 
proven to be very efficient to deal with realistic impurity 
concentrations.The average impurity spectral density 
is then given by 


Pimp(tu) — (Itn^//)avg » (5) 

n 

where (... )avg indicates the configurational average over the 
impurities. 

Figure 3 shows pimp(w) for a system with an impurity con¬ 
centration n, = 1/1800 in a cluster with 8000 impurities, dif¬ 
ferent values of the parameter V and different percentage of 
impurities on each sublattice. We consider first the less realis¬ 
tic, but simpler, case where all impurities are adsorbed on the 
Bi sites, top panels of Fig. 3. For V = 0, as in the one impu¬ 
rity case, we obtain a broad peak in pimp(w) located near the 
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FIG. 4. Average LDOS of the four non equivalent C atoms (indi¬ 
cated in the figures) in a biased BLG with impurities distributed at 
50%. The thin (red) lines are the corresponding LDOS of the pristine 
system. Upper and lower panels have V = 0.1 eV and V = -0.1 eV, 
respectively. 


renormalized energy eq. A remarkable detail is the emergence 
of a small gap for a; < 0. This gap is reminiscent of the gap 
that occurs in graphene monolayers for a finite concentration 
of impurities lying on the same sublattice.The effect is 
due to a global inversion symmetry breaking due to the differ¬ 
ent structure of the A and B sublattices. In the thermodynamic 
limit, disordered systems would not present real gaps but en¬ 
ergy windows with exponentially small DOS and it would be 
more appropriate to talk about pseudo-gaps rather than about 
real gaps. 

For a gated system with positive V = 0.1 eV the gap in¬ 
duced in the pristine BLG is partially filled by impurity states. 
Within this gap, the impurities generate a band that extends 
from the bottom of the conduction band towards the centre 
of the gap and is separated by a pseudo-gap from the valence 
band. Conversely, for V - -0.1 eV there is a narrower im¬ 
purity band close to the centre of the BLG gap separated by 
pseudo-gaps from the conduction and valence bands. These 
structures can be understood straightforwardly from the shape 
of the bound states of the one impurity case. 

The results are different if the impurities are distributed 
with the same probability on the two sublattices of the top 
layer, bottom panels of Fig. 3. For V = 0 there is no pseudo¬ 
gap on top of the valence band. For positive V an impurity 


band is formed within the BLG gap and a narrow resonance 
appears close to the bottom of the conduction band. The for¬ 
mer is due to the impurities adsorbed on the Bj sites while the 
later is due to the narrow resonance of the impurities on the 
Ai sites (see Fig. 2c). Interestingly, for negative V two sepa¬ 
rated and narrow impurity bands are formed within the BLG 
gap. Again, these bands are due to the impurities adsorbed 
on different sublattices, the lower energy one is narrower and 
comes from impurities on the Aj sublattice. 

Other impurity distributions, like the one shown at the cen¬ 
tral panel of Fig. 3 where 2/3 of the impurities are on the 
Bi sublattice, can be viewed in a first approximation as an in¬ 
terpolation between the two previous cases where the spectral 
weight of the Ai and Bi impurities change according to their 
concentration. As we show below, in gated BLG, the impurity 
bands formed within the gap of the pristine sample are bands 
of strongly localized states. 

In order to better characterize the effect of impurities on 
the electronic structure of the system we evaluate the average 
LDOS on the four non-equivalent sites of the BLG, pAi ipBi)- 
The results are presented in Fig. 4. 

Notably, in some cases, for eo > 0 and V 0, the LDOS 
of the valence band of the host BLG is almost unaffected by 
the impurities. In particular, the narrow van Hove singular¬ 
ity of the A 2 sites is essentially insensitive to the presence of 
the adsorbate. This suggests that at least in the valence band 
no strong localization effects occur with this type of impurity 
doping. As we show in next section, in gated samples, clear 
evidence of strong localization are observed for states within 
the gap and close to the impurity resonances occurring in the 
conduction band. 


B. Localization and transport properties 

To estimate the localization length ^(w) we evaluate the 
two-point correlation function |^/j(w)p, where QijioS) - 
is the retarded propagator from the impurity orbital 
at site i to the one at site j. In the localized regime this quantity 
decreases exponentially when the distance Rij between impu¬ 
rities increases."^ For large Rij (Rij > ^), the configurational 
average of its logarithm is well described by the following 
expression"^^ 

(In |^, 7 (w)|^)avg = jS - 2Rijl^(oS) - alnRij , (6) 

where a and /3 are fitting parameters. An estimation of the lo¬ 
calization length ^(w) then requires the evaluation of the im¬ 
purity propagator at large distances Rij. As discussed 

in Ref. [24], the Chebyshev polynomials method becomes 
numerically inefficient to this end. However, for long dis¬ 
tances and low energy, the propagators of the pristine BLG 
can be evaluates analytically using the continuous approxima¬ 
tion. Defining the impurity propagator matrix Q with matrix 
elements Qij{u) the Dyson equation reads 

[(w H- i0+ - eq)I - g^i. 


(7) 
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where I is the unit matrix and g is a matrix whose ele¬ 
ments are the propagators of pristine graphene, gi j(ai), be¬ 
tween C sites i and j having an impurity on top. The quantity 


tij - y^gijico) represents an effective (frequency dependent) 
hopping between impurities. The BLG retarded propagators 
take the form 


gAiiAjiR, w) = ^ cos{K ■ R) (ikiR) - h^^Kq {ik 2 R)\ 


^BZ 

A-tt 


8Ai 


i,B,j(R, W) = sm(K -R + Or) (ikiR) - (/kj/?)] , 


Q 

An 


gBu 


i,B,j(R, w) = — cos(K ■ R) (ikiR) - {ik 2 R)]. 


( 8 ) 

(9) 

( 10 ) 


Here Ky(x) is the v order modified Bessel function of the sec¬ 
ond kind, Q.BZ is the area of the first Brillouin Zone and the 
coefficients are 


{oj -v)(a) + vy 


k]iV- 


io) 


,(k2-k2) vi[k\-k^;) 


(11) 




(y H- ojf 


e. 

j 


■{kl-kl) VF{kl-kl)^ 


( 12 ) 


( (y + cu)(cu^-y^-ri) kj(V-o.) \ 

with i - 1,2, Vf is the Fermi velocity, 

kl^ = ± ^ ^I6yw + 4(w2-y2)fi. (14) 

and Qr is the polar angle of the direction of R with respect 
to the X axis chosen to be along the direction of K-K'. In 
addition, gBj^yR,oj) = gAji,B,j{-R,(^)- 

For a random distribution of impurities, we calculate the 
matrix g and obtain 0 from Eq. (7). We then take an aver¬ 
age of In \0ij{uyy for all sites i and j whose distance lies in a 
narrow window around a given value Rij. In this procedure, 
to avoid finite size effects, we take site i close to the centre on 
the cluster and neglect all sites j lying close to the edges of 
the cluster. Finally, we make a configurational average by re¬ 
peating the procedure with different impurity configurations. 
The obtained (In \0ijiu>)\^)avg versus Rij is then fitted using Eq. 
(6) to obtain the localization length ^{oj). Some of these fits 
are shown in Eig. 5 for different values of the energy at lying 
within gap of the biased BLG. In Fig. 6 the localization length 
^{a>) for the gated and ungated BLG cases and different impu¬ 
rity distributions are shown for a cluster with typical radius 
of the order of 40^,- where is the mean impurity-impurity 
distance. For the ungated system the localization length ^(w) 
presents a minimum in the conduction band for energies close 
to the renormalized energy sq of the impurity resonance. As 
the energy approaches the Dirac point from above, a> > 0, the 
localization length shows a fast increase exceeding the values 



FIG. 5. Spatial dependence of (In |^, 7 ( 6 j)P)avg inside the gap of bi¬ 
ased BLG. Top and bottom panels correspond to 100% and 50% 
of the impurities on the B[ sublattice respectively. Left and right 
columns correspond to V = 0.1 eV and V = -0.1 eV, respectively. 


for which our calculation gives reliable results (only localiza¬ 
tion lengths smaller than a fraction of the impurities cluster 
is considered). This behavior for w > 0 is qualitative similar 
to what is observed in monolayer graphene.In the valence 
band, there is a rapid increase of ^{ai) as \(x>\ increases. 

For the gated system, the impurity bands formed within the 
gap of the pristine BLG are strongly localized. In contrast, 
the states in the BLG bands tend to be much less localized, 
in particular in the valence band (consistent with the small 
sensitivity observed on the averaged LDOS, see Fig. 4). 


C. Summary and discussion 

We have analyzed the effect of diluted adatoms on the elec¬ 
tronic structure of gated and ungated bilayer graphene. The 
impurities are described as single orbital hybridized with the 
p^ orbital of one of the C atoms of the top layer. We con¬ 
sider diluted systems, typically with impurity concentrations 
n, 5 X 10 and with different distributions on the two non¬ 
equivalent sites of the top graphene layer. 

In the diluted limit studied in this work, and due to the small 
adsorption energy difference of fluorine on the two different 
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tained. In both bands the localization length shows marked en¬ 
ergy dependence with a maximum at the centre of each band. 
The states in the valence and conduction bands are much less 
localized. 

In systems with a weak energy dependence of the density of 
states and the localization length around the Fermi energy E^, 
the resistance ‘??(7’) is expected to show the Mott’s variable 
range hopping (VRH) behavior. In two dimensional system 
the VRH theory gives HiT) oc exp[(7’o/7’)5], where Tq is a 
characteristic activation temperature given by 

" k^p{E^)e{E^) ■ 


FIG. 6. Localization length in units of the average distance between 
impurities for bias voltages and different impurity distributions. Pa¬ 
rameters like in Fig. 3. Dashed lines are the average impurity spectral 
densities. 

sites, the most probable impurity distributions would corre¬ 
spond to an almost random distribution of impurities on the 
two sublattices. For the sake of concreteness, we consider the 
case of 50% of the adatoms on each sublattice, illustrated in 
the bottom panels of Fig. 6, for our following concluding re¬ 
marks. 

The first observation is that for the same impurities and the 
same concentration, the localization length is larger in ungated 
BLG than in graphene. The localization length shows a min¬ 
imum for energies close to the impurity resonance, there our 
results show that f(w) is at least two times larger in BLG than 
in graphene. The behavior of ^(«) at the Dirac point {co =» 0) 
shows a structure with a sharp minimum. In second place, the 
effect of an electric field perpendicular to the sample depends 
on the polarity of the field. The field induces a gap in the pris¬ 
tine BLG that is partially filled by strongly localized impurity 
states. However, the structure, distribution and localization 
length of these states depend on the field polarization. For 
positive V a single impurity band covers the upper part of the 
gap. There, all states are strongly localized with a maximum 
of ^(ai) at the centre of the band. The impurity spectral den¬ 
sity shows a sharp resonance at the bottom of the conduction 
band. This resonance is due to localized states that are much 
more extended than those in the gap. In the valence band the 
localization length is too large for a good estimation with the 
system size used in the calculation. For negative V two nar¬ 
row impurity bands, separated from each other and from the 
valence and conduction bands by small pseudo-gaps, are ob- 


Here a is a numerical constant (a ^ 14), p(Ep) and ^(Ep) are 
the total density of states (DOS) and the localization length 
at the Fermi energy, respectively. In biased BLG, where two 
distinct strongly localized impurities bands may exist inside 
the gap, one could expect deviations of ‘R(7’) from a single 
VRH theory. In that case, a generalization of Eq.(15) to the 
case of two narrow impurity bands might be needed. 

Finally, due to the dependence of the low energies elec¬ 
tronic structure on the polarity of the electric field, the model 
predicts a dependence of the transport properties on the sign 
of V. Such asymmetry is not clearly observed in experiments 
with fluorinated graphene.^® If all impurities where adsorbed 
on the Bi sublattice these asymmetries would be difficult to 
observe due to the similarities on the DOS and the localization 
length obtained for the two polarities, see top panels of Fig. 6. 
This scenario, however, is very unlikely. On the other hand, 
it has been reported that bilayer graphene samples grown on 
Si 02 /Si may show charge inhomogeneities with variations of 
the electronic density up to 10'*cm^. Such inhomogeneities, 
that locally shift the (electro) chemical potential in different 
regions of the sample, could also mask the asymmetries. 
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